############################################################
# Congressional Investigations and the Electoral Connection#
#           Kenneth Lowande and Justin Peck                #
# Chamber, Department, and Changepoint Analyses - 08/06/16 #
#  Please notify authors of any errors (lowande@umich.edu) #
############################################################
library(MASS)
library(MCMCpack)
library(ggplot2)

setwd("/Users/kslowande/Box Sync/Graduate School/Projects/Active/Congressional Investigations/JLEO/ms-52872-replication-files")

# Part 1: Analyses by Congress
rm(list=ls())
load("complete.Rda")
load("post1865.Rda")

# Table 3: House Models, By Congress
H1 <- glm.nb(house.x~house_div+seventeenth+h_int+cong+cong_sq+cong_cb,data=complete)
H2 <- glm.nb(house.x~house_div+seventeenth+h_int+rate_sd+exp_fd+cong+cong_sq+cong_cb,data=complete)
H3 <- glm.nb(house.x~house_div+seventeenth+h_int+rate_sd+exp_fd+war+reconstruction+cong+cong_sq+cong_cb,data=complete)
H4 <- glm.nb(house.x~house_div+seventeenth+h_int+cong+cong_sq+cong_cb,data=post1865)
H5 <- glm.nb(house.x~house_div+seventeenth+h_int+rate_sd+exp_fd+war+reconstruction+cong+cong_sq+cong_cb,data=post1865)

# Table 4: Senate Models, By Congress
S1 <- glm.nb(senate.x~senate_div+seventeenth+s_int+cong+cong_sq+cong_cb,data=complete)
S2 <- glm.nb(senate.x~senate_div+seventeenth+s_int+rate_sd+exp_fd+cong+cong_sq+cong_cb,data=complete)
S3 <- glm.nb(senate.x~senate_div+seventeenth+s_int+rate_sd+exp_fd+war+reconstruction+cong+cong_sq+cong_cb,data=complete)
S4 <- glm.nb(senate.x~senate_div+seventeenth+s_int+cong+cong_sq+cong_cb,data=post1865)
S5 <- glm.nb(senate.x~senate_div+seventeenth+s_int+rate_sd+exp_fd+war+reconstruction+cong+cong_sq+cong_cb,data=post1865)

# Part 2: Analyses by Agency
rm(list=ls())
load("complete.agency.Rda")
load("post1865.agency.Rda")

# Table 7: By Agency and Chamber Models
H1 <- glm(house~house_div*seventeenth+t_house+senate+as.factor(dp),data=complete.agency,family="binomial")
H2 <- glm(house~house_div*seventeenth+t_house+war+reconstruction+rate_sd+exp_fd+as.factor(dp),data=complete.agency,family="binomial")
S1 <- glm(senate~senate_div*seventeenth+t_senate+as.factor(dp),data=complete.agency,family="binomial")
S2 <- glm(senate~senate_div*seventeenth+t_senate+war+reconstruction+rate_sd+exp_fd+as.factor(dp),data=complete.agency,family="binomial")

# Table 8: By Agency and Chamber Models, Post-1865
H3 <- glm(house~house_div*seventeenth+t_house+as.factor(dp),data=post1865.agency,family="binomial")
H4 <- glm(house~house_div*seventeenth+t_house+war+reconstruction+rate_sd+exp_fd+as.factor(dp),data=post1865.agency,family="binomial")
S3 <- glm(senate~senate_div*seventeenth+t_senate+as.factor(dp),data=post1865.agency,family="binomial")
S4 <- glm(senate~senate_div*seventeenth+t_senate+war+reconstruction+rate_sd+exp_fd+as.factor(dp),data=post1865.agency,family="binomial")

# Part 3: Changepoint Detection
rm(list=ls())
load("complete.Rda")

# Table 10: Structural Break Models
formula <- senate.x ~ senate_div + war + reconstruction
M1 <- MCMCpoissonChange(formula,data=complete,m=1,marginal.likelihood="Chib95")
M2 <- MCMCpoissonChange(formula,data=complete,m=2,marginal.likelihood="Chib95")
M3 <- MCMCpoissonChange(formula,data=complete,m=3,marginal.likelihood="Chib95")
M4 <- MCMCpoissonChange(formula,data=complete,m=4,marginal.likelihood="Chib95")

# Table 11: Bayes Factors
print(BayesFactor(M1,M2,M3,M4))

# Figure 2: Changepoint Detection in the Senate

p0 <- data.frame(attr(M1,which="prob.state"))
names(p0)[1] <- "p0r1"
names(p0)[2] <- "p0r2"
p1 <- data.frame(attr(M2,which="prob.state"))
names(p1)[1] <- "p1r1"
names(p1)[2] <- "p1r2"
names(p1)[3] <- "p1r3"
p2 <- data.frame(attr(M3,which="prob.state"))
names(p2)[1] <- "p2r1"
names(p2)[2] <- "p2r2"
names(p2)[3] <- "p2r3"
names(p2)[4] <- "p2r4"
p3 <- data.frame(attr(M4,which="prob.state"))
names(p3)[1] <- "p3r1"
names(p3)[2] <- "p3r2"
names(p3)[3] <- "p3r4"
names(p3)[4] <- "p3r5"
names(p3)[5] <- "p3r6"

cong <- 1:80
pdata <- cbind(cong,p0,p1,p2,p3)

f1 <- ggplot() + 
    geom_vline(xintercept=64,colour="blue",linetype=2) +  
    geom_line(data = pdata, aes(x = cong, y = p0r1), colour="black") +
    geom_line(data = pdata, aes(x = cong, y = p0r2), colour="red") +
    labs(x = "", y= "") +
    theme(axis.text=element_text(size=10),
          axis.title=element_text(size=12))
f2 <- ggplot() + 
    geom_vline(xintercept=64,colour="blue",linetype=2) +  
    geom_line(data = pdata, aes(x = cong, y = p1r1), colour="black") +
    geom_line(data = pdata, aes(x = cong, y = p1r2), colour="black",linetype=2) +
    geom_line(data = pdata, aes(x = cong, y = p1r3), colour="red") +
    labs(x = "", y= "") +
    theme(axis.text=element_text(size=10),
          axis.title=element_text(size=12))
f3 <- ggplot() + 
    geom_vline(xintercept=64,colour="blue",linetype=2) +  
    geom_line(data = pdata, aes(x = cong, y = p2r1), colour="black") +
    geom_line(data = pdata, aes(x = cong, y = p2r2), colour="black",linetype=2) +
    geom_line(data = pdata, aes(x = cong, y = p2r3), colour="black",linetype=3) +
    geom_line(data = pdata, aes(x = cong, y = p2r4), colour="red") +
    labs(x = "", y= "") +
    theme(axis.text=element_text(size=10),
          axis.title=element_text(size=12))
f4 <- ggplot() + 
    geom_vline(xintercept=64,colour="blue",linetype=2) +  
    geom_line(data = pdata, aes(x = cong, y = p3r1), colour="black") +
    geom_line(data = pdata, aes(x = cong, y = p3r2), colour="black",linetype=2) +
    geom_line(data = pdata, aes(x = cong, y = p3r4), colour="black",linetype=3) +
    geom_line(data = pdata, aes(x = cong, y = p3r5), colour="black",linetype=4) +
    geom_line(data = pdata, aes(x = cong, y = p3r6), colour="red") +
    labs(x = "", y= "") +
    theme(axis.text=element_text(size=10),
          axis.title=element_text(size=12))

